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Abstract. We propose a new algorithm for generating pseudorandom 
(pseudo-generic) numbers of conformal measures of a continuous map T act- 
ing on a compact space X and for a Holder continuous potential (j> : X — > ML 
In particular, we show that this algorithm provides good approximations 
to generic points for hyperbolic rational functions of degree two and the 
potential — h\og \ T'\, where h denotes the Hausdorff dimension of the Julia 
set of T. 



1. Introduction 

Conformal measures for rational maps were introduced by Sullivan ([13]) in 
1983 following ideas of Patterson ([12]) for the case of limit sets of Fuchsian 
groups. Existence and uniqueness of such measures has been shown in [6] for 
a wide class of rational maps. These measures are in general singular with 
respect to Lebesgue measure and have no explicitly computable distribution 
function. There are a few papers dealing with the numerical computation of 
these (mostly fractal) measures (e.g. [3]), but there is no work done concerning 
the construction of generic points according to the following definition. 

Definition 1.1. Let (X, T) be a continuous dynamical system on a compact 
space X and let v be a T -invariant probability measure on the Borel field of 
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X . A point x G X is called generic if for every continuous function h G C(X) 



The convergence rate in this theorem can be arbitrarily slow, as discussed 
in [11]. It depends on the function h. A probabilistic error bound can be 
obtained from the central limit theorem or large deviation results; see jl] for a 
general description of the central limit theorem problem in dynamical systems, 
and in particular [5] about this issue for rational functions. 

Accordingly, we call a point y G X pseudo-generic for v if the equation 
(11.11) holds up to some prescribed precision or error. The construction of 
pseudorandom numbers by the linear congruential method is also based on 
the iteration theory of maps of the interval. These points well pseudo- 

generic, hence one may use the notion of pseudorandom numbers as well in 
the present situation. 

The aim of this note is to define and analyze an algorithm for computable 
pseudorandom points. In many applications, a sequence generated by the it- 
eration of a pseudo-generic point will produce points which can be viewed as 
asymptotically independent realizations of independent identically distributed 
random variables. This follows whenever the map is a weakly dependent se- 
quence of random variables. 

Let X be a compact metric space and T : X — > X be continuous. A 
conformal measure m for a continuous function <fi G C(X) is a probability 
measure satisfying 



J A 

for every measurable set A with the property that T restricted to A is invert- 
ible. This definition is equivalent to the requirement that 



for every bounded continuous function g and any measurable set A such that 
T is invertible on A and the support of g is contained in T(A). We shall call 
(11.21) the conformal equation. Examples for such measures are provided by 
rational maps (see [6] among others) or self-similar measures on fractal sets 
(see [9] among others). 

Given an invariant measure fi, equivalent to a conformal measure m, one can 
construct a pseudo-generic point by the method of least square estimation, i.e. 




n-l 



(1.1) 



k=0 




(1.2) 
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by minimizing 

E^I>Cr fe (z))- fgd>) 

gee \ k=o J / 

over a suitable subclass (5. However, the integral involved is not known and 
has to be computed by other means. It can be approximated in some cases 
using the Perron-Frobenius operator 

Ph(x)= %)exp[-<K2/)] 

T(y)=x 

and the projection to the eigenspace of the maximal eigenvalue of this op- 
erator. Here we follow another approach using a discretized version of the 
conformal equation together with a least square estimate. In this way no inte- 
gral or calculation of eigenspaces is involved in the algorithm. The algorithm 
is explained and analyzed in section 2 in general terms. 

Sections 3 and 4 are devoted to the special case of hyperbolic rational func- 
tions of degree two. We demonstrate how the algorithm is implemented in this 
case. Other cases can be investigated in a similar way. 

The algorithm requires the computation of points in X and the density of 
the invariant measure \x with respect to the conformal measure m at specific 
points. Hyperbolic rational maps on the Riemann sphere S 2 = C are charac- 
terized by the requirement that their Julia sets J(T) does not contain parabolic 
or critical points (see [I]). It is known that these maps have a unique conformal 
measure m for every Holder-continuous potential <fi : J(T) — > R (P3]), and that 
there is a unique equivalent, ergodic and T-invariant probability measure /x. 
This property will guarantee the the pseudo-generic points are approximating 
integrals with respect to the invariant measure and that the Perron-Frobenius 
operator can be used to find the density d^/dm at specific points. Unfor- 
tunately, this operator requires to calculate the Hausdorff dimension of J(T) 
which we do here numerically. Moreover, since repelling periodic points are 
dense in J(T) we are able to construct dense sets of points in the Julia set. 

It is important to remark that the algorithm needs precise numerical calcu- 
lations. We discuss this issue in Section 4. 

2. Least squares and the conformal equation 

In this section we describe an algorithm leading to explicitly computable 
pseudorandom points for a dynamical systems. We start listing the assump- 
tions we impose to hold: Let T : X — > X be a continuous map on some 
compact metric space X with metric d(-, •) and let <p G C(X) be a continuous 
function. Assume the following conditions to hold: 
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A. (1) There exists a unique conformal measure m for T and 0. 

(2) There exists a unique ergodic, T-invariant measure fi which is 
equivalent to m. 

(3) The Radon- Nikodym derivative f(x) = has a continuous ver- 
sion defined on X with modulus of continuity to : R — > M. 

B. For each n,p G N, p > n, there are finite sets X np C X, finite sets 
2l n of measurable subsets of X and continuous functions g& £ CpO 
(A G 2l„) such that 

(1) Every set A G 2l n satisfies A G er(2l n+ i) (i.e. is a union of elements 
in 2t n+ i) and every point in X lies in at most a* elements form 2l n , 
where a* is independent of n. 

(2) d :— sup^g^ max{diam(74), diam(T(74))} — > as n — > oo. 

(3) supp(^) C T(A) and < g A < 2. 

(4) The sigma fields a({g A : A G 2l n }) and cr({# A o T ■ 1 A : A G 2l n }) 
generate the Borel field of X as n — > oo. Each with a G 2l n is 
approximated arbitrarily close by a linear combination of functions 
in 2l n+ / for Z > 1 sufficiently large. 

(5) For each n, X„ iP C X UtP+ i and 

Dn iP = sup inf max d(T l (x),T\y)) — > 

xg X y&X n , v l<i<p 

as p — > oo. 

The algorithm for fixed n G N assumes the existence of m, /i, 2t n , ^ (A G 
2l n ) and the sets X njP . It proceeds as follows: 

(1) Choose za G A for A G 2l n and compute /(^a) an d f(T(zj\)). 

(2) Let p = n. For x G X njP compute /^(x) by 

£ f/(n^))E^(T^x))-/(^)E^(T fc+1 (x))e^wA . (2.1) 
agsi,-! V fc=o fc=o / 

(3) If min^gx^p < 5a*ci;(<i n ) stop and go to step 4, if not set p = 
p + 1 and continue with step 2. 

(4) Let P% = min xe x np Pn( x ) ■ Choose x* G X„ jP minimizing this expres- 
sion, i.e. 

Remark 2.1. (1) The algorithm requires to apply several subroutines explained 
by examples in the following sections: 

- Calculation of the sets X n>p and the distances D n>p for each fixed n. 

- Calculation of the sets 2l n; d n and a* . 

- Calculation of the functions gA- 

- Calculation of the density at points za andT{zA)- 
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(2) There is a simpler algorithm which may not work in general but is easier 
to implement and used in later sections. In this case all sets A have no mass 
on their boundaries (because these will be a finite union of points and the 
conformal measure has no atoms). 

The simplification puts all functions g& to indicator functions, g^ = 1t(A) 
and uses only one set X n , n . 



We need to show that the algorithm stops eventually and that the resulting 
points x* n are pseudorandom. This will be accomplished in the following two 
propositions. 



Proposition 2.2. The algorithm stops eventually. 



6 

Proof. For x G X n ^ p we have 



v 



p-i 



p-i 



< 



+ 



< 



E ( /W) E ^(T fc (x)) - f(z A ) J2 g A (T k+1 (x))e^ Tk ^ 

AeWn ^ k=Q k=Q 

E (/(TNl^A^fx))- / 
fc=o ^ 

p-i 

fc=0 

E (/(T(^))E^(T fe (a;))-/(T(z A )) / ^ 

^ fc=0 ^ 

+ /(z x ) / ^(T( M ))e^VW-/(^)E^(T fc+1 (x))e^^ 

fc=0 

E (f( T ( z A)) / ^/i- / gAfdfij 

462U \ J J J 

E (7(*a) / s(T( M ))e^/i(^) - / g A (T(u))e^f(u)fi{du 
4g2tn V Ja 

E (/(TM)E^(T fc (x))-/(T(, A )) /" ^ 
Ae2i„ ^ fc=o J 

+f{z A ) [ g A {T{u))e^)^du)- f{z A )^g A {T k +\x)) 



1/2 



1/2 



1/2 



^(r fc (x)) 



2-i 1/2 



+ 4a*u;(ci n ) 



Since /i is ergodic there exists a generic point for /x. Let z denote such a point. 
Choose x e X n , p such that d{T i {z),T i {x)) < D HjP , for each 1 < i < p. Then 



for any g A by the triangle inequality 

p-i 



I / 9Ad l i- l -Y.9A{T k {x))\ 

J P k=0 

r 1 

< | / g A dli 9A(T k (z))\ + U gA (D n ,p) 

J P i,_n 



which converges to as p — > oo, where denotes the modulus of continuity 
for the function h. Likewise 

p-i 



g A (T{u))e*%{du) - - T e^^ g A {T k+1 {x))l A {T k {x))\ 

^ fc=0 

/ g A (T(u))e^fi(du) - -^e^ Tk ^ g A (T k+1 (x))l A (T k (z))\ + u e , 9AOl 

J A P 



which tends to as well as p — > 00 (note that the support of g A oT is inside of 
A, so that the discontinuity of A is of no relevance). The indicator of A has to 
be kept here. If g A is the indicator of T(A), then g A o T{x) is one if and only 
if x is a preimage of a point in T(A). But there are more than the points in A 
as preimages and we want only those in A. 

Since 2l n is a finite set, not changing with p, we see that 

limsup —Pn — 4a*u;(<i n ), 

p— >oo P 

which implies that the algorithm stops eventually. □ 

Proposition 2.3. Let x* n and p = p{n) >n (n EN) be constructed according 
to the algorithm, and assume that A and B hold. Then for every continuous 
function g we have 

, p(n)-i 

lim — - Y, 9{T k {xD) = / gdix. 

n^ocp(n) ^ J 

Proof. Define 

p(n)-l 

^ v y fc=0 

where £ 2 denotes the point mass in z E X. Then \v n : n E N} is relatively 
compact in the weak topology of measures. Let v be an accumulation point. 
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Then v is an invariant measure and for A £ 2l n 
| J g A fdu n+l - J^ u) fg A (T(u))u n+l {du)\ 

< l/( T (^)) / 9Bdu n+l - f(z B ) j e^g B (T(u))u n+l (du)\ + o,(l) 

< /5 n+i + 0i(l). 

Therefore, letting / — > oo along a suitable subsequence, dm = /cZv satisfies 
the conformal equation, hence is conformal. Since m is unique as a conformal 
measure it equals m. Moreover, v is equivalent to m, and must be equal to /i, 
since the latter is unique as well. This shows that v n converges weakly to /i, 
proving that x* is a sequence of pseudorandom points. □ 

3. Conformal measure on a Julia set 

We describe a specific example in this section, which will be used to show 
how the general algorithm can be applied. 

Consider the rational map T : C —> C defined by T(z) = z 2 + 1. Its Julia set 
J(T) is a bounded compact set in C (under the induced (Euclidean) topology). 

A Julia set J(T) is the closure of the set of repelling periodic points and 
inverse images of T n , (n > 1) are dense as well ([I]). Thus the repelling periodic 
points and its preimages are dense inside the Julia set J(T), and every point in 
J(T) can be realized as a limit of some sequence of preimages of each repelling 
periodic point. 

The set of points y such that the iterates of y under T, T 2 , T 3 and so 
on eventually hit a fixed repelling periodic point z is dense in the Julia set. 
Therefore, it is possible to construct many points in the Julia set by taking 
preimages. This can be done in different ways. The most convenient is to 
calculate inverse branches j\ and f2 (for a quadratic polynomial) as maps 
defined on the Julia set. Then we can iterate all possible finite combinations 
/i°/2°/i°/i-j where the sequence of j\ and /2's are arbitrary choices. This 
is the naive way, since the computation creates too many data. We need for 
later purpose a certain depth of the iteration procedure, much longer than the 
forward iteration done later. In order to accomplish this one takes random 
choices of the two maps over a long string of iterations. This gives one point 
in the Julia set and one needs to estimate the errors in this calculation. 

We now discuss the discretization, i.e., a random mesh on the Julia set. 

Take a repelling periodic point z , say a fixed point: T(z ) = zq. Let us 
discretize the Julia set J(T), i.e., generate a random mesh or random lattice 
S over it, as its computational representation. 
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The inverse of T has two analytical branches, f\ and f 2 . We backward 
iterate T but select the inverse branches randomly. Let I be a large positive 
integer. Define a sample space 

Q — {1, 2}' = {uj — (co>i, (jj 2 , ■ ■ ■ , 0Ji) ■ O0i = l or 2}. 

Starting from z , a random backward iteration of I steps of T can be represented 

as 

fuJl ' ' ' fu)2 f^i 

where &x, ...,ui are chosen randomly with equal probability. Let Q be a ran- 
domly chosen set of u±, Its cardinality is the the size of the random 
mesh. Now we define a random mesh of the Julia set J(T) defined by Qq as 

S = {z : z = ■■■f Ua f ul (z ),u = (wi,w 2 , ■•• G fi }- (3-1) 

We take the Borel sets A to be small balls in J(T) centered around some 
points in S. Let us take a representative subset Sq of S and take A as balls 
centered on points in So. This is a finite family of balls and it is arranged to 
cover J{T): 

A := {A = B(z*, 5) : ball with center z* G So and radius 5 > 0}. (3.2) 

Let h denote the Hausdorff dimension of the Julia set J(T). We shall con- 
sider the conformal measure m associated to the potential hlog \T'\, which is 
a well defined Lipschitz continuous function on J(T), since T is hyperbolic. 
The transfer operator (Perron-Frobenius operator) for the the map T and the 
potential is P : C(J(T)) -> C{J{T)) defined as 

Pg{z)= £ <?(</) TO I"" (3-3) 

2/GT-i( z ) 

Iteration yields 

P n g(z)= g(y)\(T n Y(y)\- h . (3.4) 

It is known that there exists a unique conformal measure with respect to this 
potential, always denoting it m. Moreover, T is ergodic and has a unique finite 
invariant measure \x (on J(T)) that is equivalent to the conformal measure m. 

The Hausdorff dimension of the Julia set can be calculated as follows. It is 
known that 

oo 

n=0 T n (y)=x 
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converges for s < —h and diverges for s > —h. We shall use a slight variant 
of this fact to determine h: 

T«(y)=x 

converges only for s = —h. The Hausdorff dimension h is approximately 
= 1.00735 for the map z i— > z 2 + |. 

The starting point for the optimization is the defining equation for a con- 
formal measure (11.21) . We evaluate this equation for balls A C J(T). It is 
known that m has no atoms for our special quadratic map considered here, 
since the Julia set is a Jordan curve. Thus m(dB) = for open sets and we 
can construct pseudo-generic points from the equation 

m{T{A)) = [ \T'\ h dm (3.5) 

J A 

directly, where A are balls. 

Since fi ~ m, by the Radon-Nikodym theorem ^ = f(z) (density of fi w.r.t. 
m), fx(TA) = J TA fdm or 



for some point z* and z** by the intermediate value theorem since the density 
is continuous. These equations hold approximately for all z replacing z* and 
z** if the sets A and T(A) are small enough. 
Thus, on small balls A in the Julia set, we have 



f(T(z*)) n-oon k 
\ T '\ h dm = -±-r [ \T'\ h d^ 

/(**) J A 



n—l 

i— too n ^— » 



1 

— hm-EU(T fc (z))|T'(T fc (,))|\ (3.7) 



f(z*) n->oo 77, 

where z G is in a full /i— measure subset in J{T). 
Therefore, we have the fundamental equation 

1 1 nl 1 1 n ~ 1 

— —— lim -£l T(j4) (T fe (z)) = _ lim -^U(T^)) |T'(T fe (,))|*'(3.8) 
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We will find such a z approximately by the method of least squares, i.e. 
finding the minimizer for 

2 



E 



hv^ £ l nA){T\z)) - -i-i £ MT fc (,)) \T\T\ Z ))\ h 



f(T(z*))n^ v " f(z*)n 

J v v I' k=0 J y I k=0 



(3.9) 



mm 

zeS 

AeA 

where -more precisely- z* depends on the corresponding ball A. 

For calculating the density f(z*) we use the transfer operator and the well 
know equation 

/(*) = hm P n l(z) = hm V \(T n )'(y)r h (3.10) 

n—*oo n—>oo 1 ' 

yeT-™(z) 

We choose points z* G So in the discretization step and open balls around 
these points as choices of the sets A. The equation (13.101) is used to calculate 
the densities at z* and T(z*). 

We check numerically whether the minimizer in the optimization problem is 
pseudo-generic. 

For any continuous and bounded function on J{T), and for any generic point 
z in J(T), we should have 

f 1 nV 

\ gdfi = lim — } g(T h z), a.e. — [i, (3-H) 

/ ?woo n 

J k=0 

for g G L X (J(T)). In fact, in the next section, we test this for the function g = 
\z\. For z obtained in our numerical procedure, we compute lim n ^oo ^ J2k=o 9(T k (z)) 
for some large n. Repeated calculation for different sets 5", So, random choices 
and n will show that the average does not vary considerably. Choosing the 
backwards iteration randomly and So randomly may be seen as the analog of 
a seed in the construction of random numbers. 



4. Numerical results 

4.1. Determining the HausdorfF Dimension h for the Julia Set. Before 
we can hope to evaluate the Perron- Frobenius operator for a generic point in 
J(T), we must determine the appropriate Hausdorff dimension. There should 
be only one value h which allows for convergence of the limit described by 
(I3.10p to a value f(z) G (0, oo). This h will be different for each rational map 
but since we only consider one such map here, T(z) = z 2 + |, we need only 
determine the dimension of the resulting fractal (Fig. []]). 

We begin by taking the only repelling fixed point of T as a test value to 
determine h. Below is a table which shows the sequence of f n (z) values (see 
p~4j) as n — > oo. 
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0.6 
0.6 
0.4 

0.2 



-0.2 
-0.4 
-0.6 

-0.8 



Figure 1: Julia Set of the mapping T(z) = z 2 + 1/8 



Convergence of the Density Function 



z = .8356 


/i = 1.00 

fn( Z 0) fn/ fn-1 


h = 1.00735 

fn( z 0) fn/ fn-1 


h = 1.01 

fn[ Z 0) fn/ fn-1 


n = 3 
n = 5 
n = 10 
n = 15 
n = 20 


1.3029 1.0983 
1.4132 1.0299 
1.4865 1.0059 
1.5256 1.0051 
1.5644 1.0050 


1.2922 1.0935 
1.3884 1.0250 
1.4245 1.0009 
1.4258 1.0000 
1.4258 1.0000 


1.2884 1.0918 
1.3796 1.0232 
1.4028 0.9991 
1.3914 0.9982 
1.3789 0.9982 



Figure 2: Table to experimentally determine the Hausdorff dimension. 



It is clear that the approximate value h = 1.00735 allows for convergence 
of the Perron- Frobenius operator, and that values too small or too big yield 
unbounded or zero answers respectively. For the remainder of this project we 
will approximate h as such. More rigorous discussion of Hausdorff dimension 
computation can be found in [9] . If n increases the results become better in 
generally. Stratistical rigourous methods use the Grassberger and Procaccia 
correlation dimension and has been developed by Cutler and Dawson or Denker 
and Keller. The above approach is sufficient for hyperbolic maps. 

4.2. Creating the Computational Lattice. Now that we have determined 
the appropriate Hausdorff dimension, we compute the z* which define the A 
(the covering of the Julia set) and the lattice S. Given an m, the z* are all the 
points for which z = T m (z*), which can be computed directly without much 
difficulty The second block of code in Appendix A.I generates the points z* 
using zq = .8536. Logically, there are 2 m points which generate the covering 
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of the Julia Set, since there are 2 m pre-images in T m (zo). An example of a 
covering of the Julia set is below. 

Division of the Julia Set into Borel Sets 




Re(z) 

Figure 3: For m = 8 there are 256 z* points which center the circles A £ A which cover 
J(T). 



Each of the points in S is determined by an inverse iteration backwards 
from a point randomly chosen from the z*. \S\ = £, and each point is inverse 
iterated backwards N times from the initial randomly chosen point. At each 
inverse iteration only one pre-image is chosen and stored, since the other pre- 
image is of no consequence for determining the final lattice. The appropriate 
Matlab code to create the lattice is the function makelattice which is found 
in Appendix A. 2. 

4.3. Numerical Error Analysis in Random Mesh Generation. In §4.21 

we use the random mesh S to discretize the Julia set J(T). So at least, we 
would need to make sure that points in S are (approximately) inside J(T), 
when I is large enough. Numerical error here comes mainly from computing 
the inverse branches f\ = F\, f 2 = F 2 of T = z 2 + |. Let us show that this 
numerical process is stable, i.e., the total error is bounded ([TO]). 

Let F denote either one of Fi and F 2 . The error analysis for numerically 
iterating each function is similar. Let the computer's unit roundoff error be, 
for example with double precision, 10~ 16 . Let us ignore the error in computing 
the initial repelling fixed point Zq. 

Denote F(z ) be the computed value of F(z ). Then 

|F(* )-F(*o)|<|i^o)|eo, eo<10" 16 . 
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Denote Z\ = F(zq),zi = F(zq). In the following e^'s denote relative roundoff 
errors at various steps of computation. Then 

\F(~ Zl ) - F( Zl )\ = \F(z 1 )-F(z 1 ) + F(S 1 )-F(z 1 )\ 

< |F(5i) - + \F(zx) - F(z x )\ 

< \F(z 1 )\e 1 + \F / (z 1 )(z 1 -z 1 )\ 

< \F(z 1 )\e 1 + \F'(z 1 )\\F(z )\e , e x < 1CT 16 . 

As long as F and its derivative is bounded in the (bounded) Julia set, the 
right hand side of the above error estimate is upper bounded by a constant 
multiplying e < 10~ 16 . Thus when we iterate n times of F, the relative roundoff 
error ([10]) is approximately ne. The random selection between F\ and F2 does 
not change this error order. 



4.4. Efficiently Computing the Transfer Operator. Now that we have 
determined our computational lattice we must compute the transfer operator 
for obtaining f(z) which is the density as defined by the Radon-Nikodym de- 
rivative. This is a computationally sensitive segment of the procedure because 
it requires approximating a limit which becomes exponentially more expensive 
to compute. We can make this somewhat easier by noting that we consider 
only one mapping, and that T'(z) is a linear function. 

T(z) = z 2 + \ => T'(z) = 2z (4.1) 
8 

The chain rule allow us to say the following, 

(T n )'(z) = (T o T n ~ l )\z) = T'(T n - 1 (z))T'(T n - 2 (z))...T'(T(z))T'(z). (4.2) 

Below we substitute (14. ip and (14.21) into (I3.10p and use the chain rule tech- 
nique mentioned to simplify the evaluation of the transfer operator. The final 
simplification occurs because 

\uv\ = \u\\v\ Vu, «gC 

Note that the set T~ n (z) includes all 2 n (possible non-unique) values on the 
Julia set for which T n (T~ n (z)) = z. 
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/(z) = iim Yl \( Tn Y(yT h 

n—>oo ' ' 

yeT-"{z) 

= lim Yl \T\T n -\y))...T'{y)\~ h 

n—>oo 1 ' 

y£T-"(z) 

= lim V \2T n ~ 1 (y)...2y\~ h 

n. — >na * * 



lim 2 

n— >oo 



lim 2 

n— *oo 



y£T-"(z) 
hn 



E 

yeT- n (z) 



n-1 



fc=0 
n-1 



-hn 



e ni^wi 

y£T~ n (z) \k=0 



(4.3) 



We can further simplify this by noting that we are not interested in the actual 
values of T k ~ n (z) but rather only their moduli. For each value y G T~ n (z) there 
is a value — y G T~ n (z) since both the positive and negative square roots must 
be considered. Because of this we need only calculate the contribution of half 
the pre-images to the summation and then double it since \y\ = \ — y\. 

If we let T+ n denote only the positive square root pre-images of all 2 n_1 
pre-images in T~( n_1 ) (note that IT 1- ^ -1 ^ = |T_,7 n |), we can modify (14.31) 



>n-l 



Km2- hn ]II Tfc (^)l 

n— >oo ' J \ 

y£T- n (z) \k=0 



lim 2 

n— >oo 



-hn 



'n-1 



2 E \U\ Tt (y)\ 

I yeT- n (z) \k=0 

(n-l 
l\T k (y)\ 
fc=0 



(4.4) 



In the Matlab code to execute this we take advantage of the fact that we 
need only calculate and store half the moduli needed by the same logic used 
before. This is seen in graph theory and Figure H] shows how we utilize the 
fact that several branches in the tree have the same product. 

Using this logic we wrote the function dens op to approximate the limit in 
3]) to 10~ 4 accuracy. If the appropriate h value is not used, the program will 
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likely fail to converge and run in perpetuity. See Appendix A. 3 for the Matlab 
implementation. 

For those who are interested in the speed of this algorithm, ours is certainly 
not the fastest possible implementation. Each time densop is called it recalcu- 
lates pre-images which may have already been determined. In addition, new 
memory is allocated in each iteration above as well as at the start of each 
call to densop. All the code in this project is written to test the algorithms 
described above and to emphasize readability; this has resulted in a decrease 
in efficiency which will be the topic of a future project. 



4.5. Solving the Optimization Problem. Recall the optimization problem 
we need to solve: given z* and 5a (to define the Borel Set A = B(z*, 5a)) and 
n 



mm 



in > 



AeA 



f(Tz*) n 



k=0 



1 1 

f(z*) n 



n-l 



J2 1 A{T k z) \T'(T k z)\ h 



k=0 



One thing to note is that the optimization is to occur on a discrete lattice, 
S, thus we need only test a finite number of points, £, to find the solution. 
Another point of interest is that T k (z) is evaluated during construction of the 
lattice so no new function evaluations take place. Also there is no need to set 
n > £ since T £ (z) = zq for all z G S. 
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To determine if a point z is in B(z*, 5a) we simply test \z — z*\ < 5a, thus 

1 ^) = {o else ^ 

Testing if z is in _B(z*, 5t(a)) is more difficult. To do so we state that 

z G T(A), 

T- 1 ^) G A, (4.6) 
and therefore test whether either pre-image of z is in A. The result is that 

1t(A)(«)-| else (4-7) 

where T^ 1 ^) is either the positive or negative preimage of 2. 

Appendix A. 4 is the Matlab implementation of the optimization procedure. 
This code simply runs through every point in z G S and returns the point 
which minimizes (13. 9 h for summations of a given length n. It also returns the 
(3 n value described by (12. ip . 

4.6. Testing the Pseudorandom Points. In order to determine if the z 
which satisfies (I3.9P is a pseudorandom point we test its time average in equa- 
tion (13. lip . We use the simple test function g(x) = \x\ for which the integral 
can be approximated deterministically; the lhs of (13. lip is the average distance 
of points in J(T) from the origin. When all the pre-images for various values 
of m are averaged together, we see below that the result approaches the limit 
/ \z\dfj, w 1.001379. 



Integral Limit 



m 


J \z\dfi 


1 


0.853553 


5 


0.990741 


10 


1.001044 


15 


1.001369 


20 


1.001379 


21 


1.001379 



Figure 5: Table describing the approximate solution to (|3.11| . 

There are several factors which contribute to the quality of the solution and 
the complexity of the algorithm. We will assume here that we already know h 
and that m = 8 is fixed which means that the Julia set is covered by 2 8 = 256 
balls and that 5a = 2 1 ~ m is also fixed. Assuming that the z* are already 
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available (which is reasonable because the Borel sets are defined by m), the 
only parameters which affect the accuracy of the integral are 

• £ - The number of elements in S, the computational lattice 

• N - The depth of the inverse iteration on randomly chosen points z* 
to generate S 

- Recall that for z G S, T N (z) = z* for exactly one of the 2 m z*. 

• n - The point at which we truncate the least squares limit 

• a - The number of pseudorandom trajectories used to evaluate the 
integral 

We can now take a look at how changing these parameters individually 
affects the speed and accuracy of the algorithm. For all these experiments 
we have generated 10 computational lattices (ie a = 10); for each lattice 
1 pseudorandom point minimizes the least squares equation and that is the 
point whose trajectory we use to compute the integral. \x is the experimental 
mean and a is the experimental standard deviation 

We can see from Figure [6c] that increasing N has the effect of decreasing 
the standard deviation of the estimator. Figure [7] is a graphic depiction of 
this. It appears in Figure [6b] and Figure [6a] that increasing n and £ without 
changing N causes no improvement in a. This leads us to believe that the 
driving force behind accuracy is N for which the complexity of the algorithm 
increases linearly. 

There are limitations to this algorithm because it requires storage of N+l x£ 
terms: this is done to prevent the loss of accuracy from N applications of 
T on the elements of S. Unfortunately, since there is only one point in S 
which minimizes (13.91) there is only one pseudorandom trajectory chosen per 
S generated. 

One possible future improvement to this algorithm may be to use several tra- 
jectories whose value in (13. 9p are close to optimal but not the exact minimum. 
Their contribution to the estimator can be weighted according to their distance 
from the optimal value. This would allow the use of multiple trajectories from 
the same S and not require a versions of S for a trajectories. 
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(a) Fixed N = 16000, n = 


100. 


(b) Fixed N = 


16000, 1 = 
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Figure 6: The effect of varying N, n and 



Appendix A - Matlab Algorithms 

The first file is a script which calls the other functions to generate pseudorandom numbers 
on the Julia Set for the mapping T(z) = z 2 + 1/8. Here is a list of important parameters: 

• zO - A repelling periodic point which is the start of the inverse iterations 

• m - 2 m Borel sets are used to cover the Julia Set 

• zstar - The centers of the Borel sets. These are the 2"m pre-images in the set 
T~{-m}(z ) 

• S - The discretization of the Julia Set 

• ell - The number of points in S 

• N - All points in S are pre-images in the set T"{— (m + N)}(zq) 

• n - Summations in the optimization equation are of length n 

• alpha - The required number of pseudorandom points 

• h - The Haussdorff dimension 



A.l MAINSCRIPT.M 
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Convergence of Confidence Interval 
1 .005 1 ■ ■ ■ — ■ ■ 




Figure 7: Solid - ji ± cr, Dashed - /i, Dotted - true solution 



% Here we consider the mapping T( z) = z "2+1/8 
T=inline ( 'x.*x + .125 ' ) ; 

z0=fsolve (@(z) T(z)-z,.8); % approximately zO =0.85355339203135 
m=8; 
ell =100; 
N=32000; 
n=100; 
alpha =30; 
h = l. 00735; 

indmin = zeros ( alpha , 1 ) ; 
beta_n = zeros ( alpha , 1 ) ; 
caverage = zeros ( alpha , 1 ) ; 

% This line makes the random number generator start with the 
% same seed always. The line below will randomize the seed. 
rand( 'state ' ,0) ; 

% rand ( ' state ' , sum (100* clock )) ; 

% We use inverse iteration to create the lattice. We must 
% determine the points in zstar and define the balls which 
% cover the Julia Set . 
zstar =z0*ones(2"m,l); 



21 



for bdec = l:2"m 

bcode=dec2bin ( bdec — 1 ,m) ; 
for i=l:m 

zst ar ( bdec ) = ( — 1) * ( bcode ( i)== ' 1 ' ) * sqrt (zstar (bdec) — .125); 
end 
end 

% The function densop finds the density at zstar . 
% This does not have a time limit , so if the script is 
% hanging , densop is a likely source . Both the density of 
% the points in zstar and their images T( zstar) are found. 
fzstar=zeros (2~m, 1 ) ; 
fTzstar=zeros (2~m, 1 ) ; 
for k = l:2:2"m 

fzstar (k) = densop(zstar (k) ,h); 

f z st ar (k+l)= densop (zstar ( k + 1) ,h ) ; 

fTzstar(k)=densop(T( zstar (k ) ) ,h); 

fTzstar (k+l)=fT zstar (k) ; 
end 

% Now use makelattice to form the discrete Julia Set . The 
% seeds for the inverse iteration of points on the lattice 
% are randomly chosen from zstar. The computational lattice 
% is S(: ,N) . The final column is an extra inverse iteration 
% for evaluating the fundamental equation. It is important 
% to note that the size of S is [ ell ,N+1] not [ell,N]. 

% After that we use the optimization function which will return 
% the solution to the least squares problem detailed earlier. 
% The function below returns imin , the index of the point in S 
% which is the solution , and minival which is the residual . 
% We should have beta_n /n"2<2*dn * | A_n j * LipConst as described 
% in the earlier paper, 
for i=l:alpha 

S=makelattice ( ell ,N,m, zstar); 

[imin(i),beta_n(i)] = opteval(n,S, zstar , fzstar , fTzstar ,m, h); 
end 

% Now we test the ensemble averaging . Any LI function can be 
% used to test the pseudorandomness of the points found by the 
% optimization . The function must be able to accept vector 
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% arguments, ie using .* instead of just * for multiplication. 
g= inline ( 'abs(x) ' ); 
for i=l:alpha 

caver age ( i )=mean( g ( S ( imin ( i ) , 1 :N) ) ) ; 
end 

A. 2 MAKELATTICE.M 

function S=makelattice ( ell ,N,m, zstar) 
% function S=makelattice ( ell ,N,m, zstar) 

% This function makes the computational lattice that represents 

% the Julia Set . The rest of the matrix values are 

% the trajectories taken backwards from a random selection of 

% points on the zstar grid. S(:,N) is the computational 

% lattice , S(:,l) is ell randomly chosen points from zstar. 

% S(: ,N+1) is a preimage of S(: ,N) which is needed for 

% optimization to test if S(: ,N) is in TA. 

S = 2*(rand(ell ,N+1)>.5)-1; 

S ( : , 1 ) = zstar(ceil (2"m*rand( ell , 1 ) ) ) ; 

for j=2:N+l 

S(: ,j) = S(: ,j ).*sqrt(S(: , j -1) -.125); 

end 

A. 3 DENSOP.M 

function f z=densop ( zO , h) 

% function f z=densop ( zO , h) 

% This function computes the density for the 

% mapping T( z)=z "2 + 1/8. It does this using a limiting sequence, 
fval =[0 , -1]; 
c = l; 
cc = 1 ; 
pq(l) = z0; 
pn(l) = abs(pq(l)) ; 
k = 2; 

while abs(fval(2)-fval(l))>le-4 
c = [c,2"(k-l)]; 
cc=cumsum( c ) ; 
pq=[pq, zeros (1 ,c(k) )] ; 
pn=[pn , zeros ( 1 , c (k ) ) ] ; 
fval(l)=fval (2); 
fval (2) = 0; 
for j=c(k) :2: cc(k) 
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pq(j)=-rj*sqrt(pq(fix(j/2))-.125); 

pq( j+i)=-pq( j ) ; 

pn(j)=abs(pq(j )); 
Pn( j+l)=pn( j ) ; 
end 

kk=c(k) ; 
while kk<cc ( k) 
kj=kk; 

while kj ( length ( kj )) >2 

kj=[kj , kj ( length (kj))/2-(mod(kj( length (kj)) ,4) >0)]; 
end 

fval (2)= fval (2) + (prod (pn( kj ))) * -h ; 
kk=kk + 2; 
end 

fval (2)= fval (2)*2 A (-h*(k-l) + l); 
k=k + l; 
end 

fz = fval (2) ; 

A. 4 OPTEVAL.M 

function [ indmin ,minival] = opteval(n,S,zstar , f zst ar ,fTzstar,m, h) 
% function [ indmin ,minival] = opteval(n,S,zstar ,fzstar , fTzst ar ,m, h ) 
% This needs the computational lattice , the transfer operator 
% evaluated on the lattice , and the size of the Borel Set A 
% around zt . n is the limit truncation which can not be greater 
% than N, and h is the Hausdorff dimension. 

% All this function does is run through the lattice and calculate 
% the optimization equation at each point . It returns the lowest 
% value. S is the group of trajectories which yield the 
% computational lattice. Specifically S(:,N) is the lattice, 
[ell ,N] = size(S); 

N=N-1; % Recall size (S) = [ ell ,N+1] although S(:,N) is the lattice. 
% The radius of the balls which cover the region is related to m 
delta_A=2~(-m+l); 

% The lhs part of the summation will find whether either of the 
% preimages of S are in A. This is equal to asking if S is in 
% T(A) . The rhs part of the summation tests whether S is in A, 
% and then adds the appropriate values. The rest is just 
% evaluating the optimization equation. 
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Asum = zeros (ell ,1); 
for j = l:2~m 

spS = sparse (abs (S (: ,N—n :N)— zstar ( j ))< delta_A ) ; 

rhs = 2*sum((spS.*abs(S(: ,N-n :N) ) ) . ~ h , 2 ) ; 

lhs = sum((S(: ,N+l-n:N+l)-zstar ( j)<delta_A) + ... 

(-S(: ,N+l-n:N+l)-zstar (j)<delta_A) ,2); 
Asum = Asum+( lhs / fTzstar (j) — rhs / fzstar (j ) ) . A 2; 

end 

[ minival , indmin ] = min(Asum/N) ; 
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